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A fluid in equilibrium in a finite volume V with particle number A at a density p = N/V exceeding 
the onset density p/ of freezing may exhibit phase coexistence between a crystalline nucleus and 
surrounding fluid. Using a method suitable for the estimation of the chemical potential of dense 
fluids we obtain the excess free energy due to the surface of the crystalline nucleus. There is neither 
a need to precisely locate the interface nor to compute the (anisotropic) interfacial tension. As a 
test case, a soft version of the Asakura-Oosawa model for colloid polymer-mixtures is treated. While 
our analysis is appropriate for crystal nuclei of arbitrary shape, we find the nucleation barrier to be 
compatible with a spherical shape, and consistent with classical nucleation theory. 


Nucleation of crystals from fluid phases and their 
subsequent growth is one of the most important phase 
transformations in nature CHS]; applications range from 
ice crystal formation in the atmosphere, to metallurgy, 
nanomaterials, protein crystallization, etc. Despite its 
overwhelming importance, crystal nucleation still is only 
poorly understood. 

For the nucleation of a liquid drop from supersatu¬ 
rated vapor, clearly the average nucleus shape is spheri¬ 
cal. Only the curvature dependence of the interfacial ten¬ 
sion Hi] presents a stumbling block for the prediction 
of nucleation barriers. Unlike interfaces between fluid 
phases, the crystal-fluid interface tension 7 (n) depends 
on the orientation of the interface normal n relative to the 
crystal lattice axes [inUT^ . For isotropic 7 the nucleus is 
a sphere of radius R (volume V = 47 ri?^/ 3 ) and its sur¬ 
face excess free energy is Fgurf = 4-71 R^j = 
with Aiso = (SGtt)^/^. For crystals the term AisoJ is 
replaced by a complicated expression, 

Fsurf(V) = I j{n)dsV^/^ = . (1) 

Tw 

Here Ayy is the surface area of a unit volume whose shape 
is derivable from 7 ( 7 ?) via the Wulff construction [TDl - 
[ 12 ] , and the average interface tension 7 is defined as 
7 = S lin)ds. 

In the classical nucleation theory m, the formation 
free energy of a nucleus is written in terms of volume and 
surface terms as 


AF = -{p,-pi)V + Fsurf{V) . (2) 

Here Pc is the pressure in the crystal nucleus and pi in 
the (metastable) liquid phase surrounding it. In the ther¬ 
modynamic limit, the configuration with one nucleus on 
top of the free energy barrier in the metastable phase is 
a saddle point in configuration space. The condition for 
(unstable) equilibrium, d{AF)/dV = 0, then yields the 





FIG. 1. Schematic plot of the chemical potential p. versus 
density p for a system undergoing a liquid-solid transition in 
a finite box volume Vbox with periodic boundary conditions. 
(A plot of pressure p vs. p would qualitatively look just the 
same). Due to interfacial effects, non-negligible in finite sys¬ 
tems, the isotherm deviates from p = Pcoex in the two-phase 
coexistence region, pf < p < pm- The features in the curve 
(kinks in reality are rounded due to fluctuations) are due to 
transitions between the different states shown in the figure 
via snapshots of the simulated generalized Asakura-Oosawa 
model. Only the part where the solid phase is the minority 
phase is discussed. For further explanations cf. text. 
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Even if V* is large enough so that correction terms to 
Eq.i can be neglected, the application of Eq. is diffi¬ 
cult due to lack of knowledge on Aw and 7 . This lack 
of knowledge has hampered the comparison of observed 
nucleation rates [T5HTB] (and the barriers extracted from 
them) and simulations [T7IED] where AF* was estimated 
directly by biased sampling methods. These comparisons 
were made for suspensions of (hard sphere-like) colloidal 
particles; the large size of the colloids has the advantage 
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FIG. 2. Normalized pressure p — pa^/ksT plotted vs. pack¬ 
ing fraction rj = piva^fQ of the colloids, for the effective AO 
model (henceforth denoted as Eff AO, asterisks) and its soft 
version (soft Eff AO, squares). Curves are guide to the eye 
only. These data were obtained from simulations of homoge¬ 
neous liquid and solid (fee) phases, while the pressure where 
two-phase coexistence occurs was found from the “interface 
velocity method” |37| . namely p = 8.44 ± 0.04 (soft Eff AO) 
and p = 8.06 ± 0.06 (Eff AO). The coexistence packing frac¬ 
tions are 77 / = 0.495 (1) and pm ~ 0.636(1) for the soft Eff 
AO case. The insert compares the potentials of the Eff AO 
(which is singular at r — ac = 1) and soft Eff AO models. 


to allow direct microscopic observations of crystal-liquid 
interfaces [2T] and nucleation events [22l |23] . Since ki¬ 
netic processes for colloids are many orders of magnitude 
slower than for small molecules, colloids are model sys¬ 
tems for the study of the liquid-solid transition [531 US] j 
and well suited to separate nucleation from the subse¬ 
quent crystal growth. 

However, to elucidate the persisting discrepancies be¬ 
tween simulations and experiments one needs to know 
more about the theoretical nucleation barriers: How large 
must V* be so that Eq.j^is a good approximation? What 
is the physical origin of corrections to AF* {Eq. H and 
their magnitude? Is it legitimate to assume a spherical 
shape of the nucleus, in spite of its crystalline structure? 
And so on. Understanding the general conditions under 
which the classical description {Eqs. in holds will be 
useful to understand liquid-solid transitions in condensed 
matter in general. 

In the present letter, we address these issues, and show 
how both V* and AF* can be obtained, considering the 
equilibrium of the system at fixed finite particle number 

in a finite simulation box Vbox- For a suitable range 
of density p = N/Vbox, the equilibrium between the crys¬ 
talline nucleus and surrounding fluid is perfectly stable. 
We explain how both V* and Pc — Pi can be estimated 
directly and accurately. Using then AF* = {pc—pi)V*/2 
{Eq. [^, the need of dealing with 7 (h) and use of Eq. 
is bypassed. So we do not need to assume anything on 
the shape of the nucleus. 

Thus, the central idea of the present work is to explore 
the deviations from phase coexistence in the thermody¬ 


namic limit (where the chemical potential p, = Pcoex and 
the pressure p = Pcoex for all densities from the onset 
density of freezing pj to the onset density of melting 
caused by finite size. Thus, the part of the isotherm in 
Figj^ corresponding to the homogeneous fluid for finite 
volume Vbox exceeds p / and continues up to the “droplet 
evaporation condensation transition” [26] at pi , where for 
the first time a crystalline droplet in the system becomes 
stable. Note that this transition is a sharp phenomenon 
only when Vbox —>■ 00 (and then pi —>■ pf, consistent 
with the lever rule m)- At a second special density p 2 
the “droplet” changes its shape from compact to cylin¬ 
drical (stabilized by the periodic boundary conditions). 
At about p = p 3 a slab configuration, separated from the 
fluid by two planar interfaces, appears (Fig.l). In this re¬ 
gion p = Pcoex and p = Pcoex holds true also in the finite 
system, if the linear dimensions in the directions parallel 
to the planar interfaces are chosen such that the crystal 
(at density pm) is commensurate without any distortion. 
The analogous behavior for vapor to liquid transitions is 
well studied jO] 1281130] . Here we show that the descend¬ 
ing part of the p(p) and p(p) isotherms can be used to 
extract information on Fsurf, V* and AF* for the liquid- 
solid transition as well. 

In the snapshots the particles in the fluid region are 
shown in blue, in the crystal in red color, using the av¬ 
eraged Steinhardt local bond order parameters 
to distinguish the character of the phases (see B ef. [05] 
for definitions and implementation details). Particles in 
the interfacial region, for which this classiheation yielded 
ambiguous results, are shown in green color. The face- 
centered cubic (fee) packing of the crystal is clearly seen, 
and the cross section through the “droplet” also suggests 
that the shape may non spherical. 

The model of our simulations qualitatively describes 
colloid polymer mixtures [33ll36] . In the Asakura- 
Oosawa (AO) model [33], colloids are described by hard 
spheres of diameter CTc, polymers as soft spheres (which 
may overlap each other without energy cost) of diameter 
CTp. Of course, the mutual overlap of colloids and poly¬ 
mers is also strictly forbidden. Polymers create the (en- 
tropic) depletion attraction between colloids [53] ; varying 
the size ratio q — Opjac and the polymer density one can 
tune the phase diagram [341436] and interfacial proper¬ 
ties [371135]. A useful feature of this model occurs for 
q < q* = 0.154 [33[[3n]: then one can integrate out the 
polymer degrees of freedom exactly, and one is left with 
an effective pairwise potential, which is attractive in the 
range ac < r < ac + Up (and zero for r > ac + Cp), but 
infinitely repulsive for r < ac- The strength of the po¬ 
tential of this “effective” AO model is controlled by the 
fugacity Zp of the polymers [39] (Fig. 2, insert). 

However, it is computationally more convenient to re¬ 
place the Eff AO model by a similar but continuous po¬ 
tential, the soft Eff AO model [39] (Fig. 2, insert). For 
this model the pressure (in the fluid phase) is straight¬ 
forwardly obtained in the simulation from the Virial ex¬ 
pression [SailO], while for the Eff AO model due to the 
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discontinuity at r = CTc this is very cumbersome |38j . 
Fig. 2 shows that the variation of p with 77 is very simi¬ 
lar for both potentials. Since real colloids never are de¬ 
scribed by hard spheres precisely m, nor are polymers 
precisely modeled by ideal soft spheres [ 12 ], a quantita¬ 
tively accurate modeling of real systems anyway cannot 
be attempted. The soft Eff AO model is proposed here 
as a coarse-grained qualitative model of colloid-polymer 
mixtures which is practically useful in a simulation con¬ 
text. 

Using the Virial expression the pressure pi of the liq¬ 
uid in the region surrounding the crystal nucleus in Fig. 1 
(far away from the interfacial region) can be readily mea¬ 
sured, but obtaining pc inside the nucleus for small nuclei 
is not reliably possible. It is necessary to base the anal¬ 
ysis of the two-phase equilibrium in Vbox on the chem¬ 
ical potential p, because p is strictly constant in equi¬ 
librium also in a spatially inhomogeneous situation. But 
the standard particle insertion method gniiig does not 
work at high packing fractions pc near pm- Thus, we have 
extended an approach [44] to sample the chemical poten¬ 
tial of a dense fluid by studying a system where walls 
are present; using a soft wall that reduces the density 
suitably such that there particle insertion works (Fig. 3). 
Of course, it is important to choose large enough so 
that outside of the range of z, for which the walls affect 
the density prohle, actually a constant density is reached. 
Fig. 3 demonstrates that in this way the chemical poten¬ 
tial can be obtained accurately even for 77 > pf. The pres¬ 
sure p (computed in the region where p(z) = Pbuik = con¬ 
stant) agrees with the corresponding bulk data of Fig. 2. 

Now we exploit the fact that p is constant through¬ 
out the system also when a crystalline nucleus is present 
(Fig. 1): the chemical potential in the fluid Pf{pi) equals 
that of the crystal nucleus pdpc)- From pdPcoex) = 
PiiPcoex) = Pcoex we readily find, using the expansions 

TT 1 

/^c(Pc) ~ l-^coex 77 (.Pc Pcoex) •> (4) 

6 r]m 

TT 1 

Pl(Pl) ~ Pcoex H” "T (pi Pcoex) i (5) 

6 Pf 

that {pc - Pcoex)pf = {pi - Pcoex)pm- Since we have 
recorded both functions pi{p) and pi{p), we also know 
Pi{pi) and hence can verify that the data indeed fall 
in the regime where the linear expansion, Eq. [^ holds. 
Finding pdpd) via thermodynamic integration (using 
pdPcoex) = PliPcoex) as Starting point), we have verified 
that Eq. @ also introduces only negligible errors. 

The two-phase equilibrium of a crystalline droplet sur¬ 
rounded by fluid has been studied for three system sizes, 
keeping the number of colloids in the simulation box fixed 
(at N = 6000, 8000, and 10 000, respectively) and vary¬ 
ing Vbox and hence p = N/Vbox- In thermal equilibrium, 
we then have a finite-size variant of the lever rule 

pVbox = Pi{pi){Vbox - V*) + pdPc)V* (6) 

While for Vbox —>■ oo we would have Pi = Pc = Pcoex 
and pdPcoex) = Pm, in the finite system pi,Pc and the 
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FIG. 3. (a) Illustration of the method to compute the chemi¬ 

cal potential of a very dense fluid, using a. LxLxLz slab geom¬ 
etry, with a soft wall at z = 0 and a hard wall at 2 = L 2 = 30 
(lengths being measured in units of CTc, L = 7, and 4 choices 
of N are used, N = 750, 950,1100 and 1250, respectively). In¬ 
sert shows p in units of ksT as a function of 2 , for the 4 choices 
shown, over the regions of 2 where particle insertion works, 
(b) Chemical potential p (in units of ksT) plotted vs. 77 , for 
different choices of L and Lz, as indicated, to show that finite 
size effects are negligible. The data labeled by “Widom” at 
not so large p are obtained by the standard particle insertion 
method for homogeneous bulk systems. Arrows on abscissa 
and ordinate indicate 77 / and pdVcoex)/ksT, respectively. 


corresponding packing fractions differ from their coex¬ 
istence values. Initializing the simulation by putting a 
crystal of about the right volume V* and about the right 
choice for pdPc) in the box, after a long period of equi¬ 
libration we measure both pi and pi{pi) in the fluid re¬ 
gion (far away from the crystal) and verify (from the 
data of the bulk simulation. Fig. 3b) that equilibrium has 
been reached. Since we know also the chemical potential 
(Mc(Pc) = Piipi) is constant), we can obtain pc and also 
Vc(Pc) and hence Eq. determines V* unambiguously. 

Fig. 4 shows the data for Ap — pi — pcoex versus p. 
Actually when we use the chemical potential pi{pi) from 
Fig. 3b and obtain Ap from Eq. [^ the data are precisely 
reproduced, which just is a consistency check. From sim¬ 
ulations determining 7 for interfaces parallel to 111,110 
and 100 planes [IS] it is found that 7 ( 71 ) depends only 



















4 


very weakly on ft. For comparison with classical nucle- 
ation theory, we neglect the dependence on ft and take 
7 « 7 m « 7 (n) = 1.013 [3S]. Assuming a spherical 
shape V* = 47ri?*^/3 we find Ap = {2j/R*)/{r]m/'r]f~^)- 
Using the observed values of V* one then obtains a 
prediction for the curves Ap{r]). We find that these 
predicted curves fall slightly below the actual observed 
data. They can be brought in good agreement if they are 
rescaled by a constant factor of c = 1.07. This small en¬ 
hancement can be due to the ratio A/Aiso or errors in the 
estimation of "/(ft). Unexpectedly, we hence find that for 
our model of colloid polymer mixtures the assumption of 
a spherical nucleus shape works rather well, but it would 
not be needed to predict the nucleation barrier. Using 
Eq-d knowledge of Pc — Pi and V* suffices to predict 
AF*. One can expect, however, that significant deriva¬ 
tions from spherical nucleus shape will appear for large 
rjp in our model, where the fluid is a vapor-like phase, 
and 7 (n) will depend more strongly on ft. Gratifyingly, 
Fig. 4b shows that the three choices for N superimpose 
to a common curve, so in the shown regime finite size 
effects essentially are negligible. 

In summary, we have shown that for the liquid-solid- 
transition a description of nucleation barriers in terms 
of the classical nucleation theory holds, at variance with 
studies of nucleation with hard sphere-like colloids [ni- 
unmii. However, we feel the latter studies are inconclu¬ 
sive, do to their use of too large pi (0.53 < pi < 0.57) 
where the slowing down due to the kinetic prefactor of 
the nucleation rate matters m- While the range of 
AF* in Fig. 4b corresponds to pi/pf — 1 < 0.06 the 
range of the experiments in Fig. 4b would correspond 
to 5 < AF* < 10 only. 

Analyzing finite size effects on phase coexistence, both 
V*, pi, Pc and the chemical potential for this stable 
two-phase coexistence in a finite simulation box can be 
reliably estimated. The numerical results also clearly 
show that in the regime where AF* > 80 the relation 
AF* oc y*’^!'^ holds precisely, as visible from the fit 
in Fig.4b; thus we have verified that classical theory of 
homogeneous nucleation for crystals is accurate, in this 
regime of barriers, provided one takes into account that 
the nucleus shape is in general nonspherical. However, 
since the two straight lines in Fig.4b almost coincide, the 
spherical approximation is shown here to be almost per¬ 
fect. Since crystal faces in contact with a dense fluid are 
frequently atomically rough, the spherical approximation 
is expected to be quite good generally, in particular for 
somewhat smaller nuclei, for which the nucleation rates 
also would be larger. 
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FIG. 4. (a) Pressure difference Ap = pi — Pcosx between the 

pressure pi in a fluid surrounding a crystal nucleus of finite 
size and the coexistence pressure, plotted versus the average 
packing fraction p in the simulation box, for particle number 
N = 6000, 8000 and 10000 (symbols, from bottom to top). 
Curves show the formula Ap — {2pjR*)cl{pm/pf — 1) with 
c = 1.07, extracting R* from the assumption of a spherical 
nucleus {V* = A-kR*^ /Z) and taking 7 « 7111 « 1.013 |4^ . 
(b) AF* computed fromp;,pc and V* (using Eqs.[^|^and|^ 
plotted vs. y*'^A^ straight line is Eq.(3) with A — Aiso and 
7 = 7 iii- The broken line is a fit illustrating AF* oc 
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